# load packages
library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.1 ✔ purrr 0.3.2
## ✔ tibble 2.1.3 ✔ dplyr 0.8.3
## ✔ tidyr 1.0.0 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0
## ── Conflicts ───────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(sf)
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(leaflet)
library(rcartocolor)
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:tidyr':
##
## extract
library(dplyr)
library(tmap)
library(rasterVis)
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
##
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
##
## layer
library(maptools)
## Checking rgeos availability: TRUE
library(fasterize) # fasterize is high-performance replacement for the rasterize
##
## Attaching package: 'fasterize'
## The following object is masked from 'package:graphics':
##
## plot
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(rcartocolor)
library(rmapshaper)
## Registered S3 method overwritten by 'crul':
## method from
## as.character.form_file httr
## Registered S3 method overwritten by 'geojsonlint':
## method from
## print.location dplyr
library(stringr) # for working with strings (pattern matching)
library(tidyr) # for unite() and separate()
library(spData)
library(velox) # faster extraction in raster
library(hrbrthemes)
detach("package:dplyr", unload = TRUE)
## Warning: 'dplyr' namespace cannot be unloaded:
## namespace 'dplyr' is imported by 'broom', 'tidyr' so cannot be unloaded
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:gridExtra':
##
## combine
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
You should understand and have control over the geometry column in sf objects and the geographic location of pixels represented in raster
The sf pacakge provides st_simplify(), which uses the GEOS implementation of the Douglas-Peucker algorithm to reduce the vertex count. st_simplify() uses the dTolerance to control the level of generalization in map units
par(mfrow = c(1, 2))
seine_simp <- st_simplify(seine, dTolerance = 2000)
plot(seine)
plot(seine_simp) # 2000 m
the resulting is a copy of oringal seine but with few vertices. This is, apparent, with the result being visually simply and cosuming less memeory than the original objects
map_dbl(list(seine, seine_simp), object.size)
## [1] 17304 8320
GEOS assume that the data is in a projected CRS and this could lead to unexpected results when using a geographic CRS.
us_states_projected <- st_transform(us_states, 2163)
# st_simplify works equally well with projected polygons
us_states_simp1 <- st_simplify(us_states_projected, dTolerance = 100000) # 100km
plot(us_states_simp1 %>% st_geometry())
limitation with st_simplify is that it simplifies objects on a per-geometry basis. This means the topology is lost, resulting in overlapping and holy areal units. ms_simplify from rmapshaper provide an alternative that overcome this issue.
# proportion of points to retain (0-1, defafault .05)
us_states_projected$AREA <- as.numeric(us_states_projected$AREA)
us_states_simp2 <- rmapshaper::ms_simplify(us_states_projected, keep = 0.01, keep_shapes = T)
plot(st_geometry(us_states_simp2))
Centroid operations identify the center of geographic objects. Like statistical measures of central tendency (including mean and median definations of average), there are many wasy of define the geographic center of an object. And sometimes the geographic centroid falls outside the boundaries of the aprent objects. In such case points on surface operations can be used to guarantee the point will be in the parent object. In this case use st_point_onsurface()
nz_centroid = st_centroid(nz)
## Warning in st_centroid.sf(nz): st_centroid assumes attributes are constant
## over geometries of x
nz_pos <- st_point_on_surface(nz)
## Warning in st_point_on_surface.sf(nz): st_point_on_surface assumes
## attributes are constant over geometries of x
seine_centroid <- st_centroid(seine)
## Warning in st_centroid.sf(seine): st_centroid assumes attributes are
## constant over geometries of x
seine_po <- st_point_on_surface(seine)
## Warning in st_point_on_surface.sf(seine): st_point_on_surface assumes
## attributes are constant over geometries of x
plot(st_geometry(nz))
plot(st_geometry(nz_centroid), col = "black", add = T)
plot(st_geometry(nz_pos), col = "red", add = T)
plot(st_geometry(seine))
plot(st_geometry(seine_centroid), col = "black", add = T)
plot(st_geometry(seine_po), col = "red", add = T)
BUffers are polygons representing the area within a given distance of a geometric feature: regarding of whether the input is a point, line or polygon. Unlike simplifications (which is often used for visualization and reducing file size) buffering tends to be used for geographic data analysis (e.g.: calculate does different race people live close to school district or not)
st_buffer requireds at least two arguments: an input geometry and a dstance
seine_buff_5km <- st_buffer(seine, dist = 5000)
plot(seine_buff_5km)
Affine transformation is any transformation taht preserves lines and parallelism. However, angles or lengths are not necessarily preserved. Affine transformation inlcudes, shfiting, scaling and rotation
Rotation matrix
$ = $
nz_sfc <- st_geometry(nz)
# shifts all y-coordinates by 100,000 meters to the north,
#but leaves the x-coordiantes untouched
nz_shift = nz_sfc + c(0, 100000)
plot(nz_sfc)
plot(nz_shift, col = "red", add = T)
# scaling enlarge or shrinks objects by a factor, be applied either globall or locally
nz_centroid_sfc <- st_centroid(nz_sfc)
nz_scale = (nz_sfc - nz_centroid_sfc) * 0.5 + nz_centroid_sfc
plot(nz_sfc)
plot(nz_scale, col = "red", add = T)
# rotation
rotation <- function(a) {
r = a * pi / 180 # degrees to radians
matrix(c(cos(r), sin(r), -sin(r), cos(r)), nrow = 2, ncol = 2)
}
nz_rotate <- (nz_sfc - nz_centroid_sfc) * rotation(30) + nz_centroid_sfc
plot(nz_sfc)
plot(nz_rotate, col = "red", add = T)
spatial cipping is a form of spatial subsetting that involves chagnes to the geometry column of at least some of the affected features. Clipping can only apply to featrues that are more complex than points: lines, polygons and their “multiple” equivalents
b <-st_sfc(st_point(c(0, 1)), st_point(c(1, 1)))
b <- st_buffer(b, dist = 1)
plot(b)
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"))
if ypu want to select not one circle or the others, but the space covered by borth x and y and rest of “Venn” diagram representing x and y
x <- b[1]
y <- b[2]
x_and_y <- st_intersection(x, y)
plot(b)
plot(x_and_y, col = "lightgrey", add = T)
# in y but not x
plot(b)
plot(st_difference(y, x), col = "lightgrey", add = T)
text(x = 0.5, y = 1, labels = "st_difference(y, x)")
# both x and y
plot(b)
plot(st_union(y, x), col = "lightgrey", add = T)
text(x = 0.5, y = 1, labels = "st_union(y, x)")
# not both in x and y
plot(b)
plot(st_sym_difference(y, x), col = "lightgrey", add = T)
text(x = 0.5, y = 1, labels = "st_sym_difference(y, x)")
box <- st_as_sfc(st_bbox(st_union(x, y)))
set.seed(2017)
p <- st_sample(x = box, size = 10)
plot(box)
plot(x, add = T)
plot(y, add = T)
plot(p, add = T)
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"))
the logical operator way would find the points inside both x and y using a spatial predicate such as st_intersects
sel_p_xy <- st_intersects(p, x, sparse = F)[, 1] &
st_intersects(p, y, sparse = F)[, 1]
p_xy1 <- p[sel_p_xy]
p_xy2 <- p[x_and_y]
identical(p_xy1, p_xy2)
## [1] TRUE
Spatial aggregation can silently dissolve the geometries of touching poloygons in the same group.
regions <- aggregate(x = us_states[, "total_pop_15"], by = list(us_states$REGION),
FUN = sum, na.rm = T)
regions2 <- us_states %>%
group_by(REGION) %>%
summarize(pop = sum(total_pop_15, na.rm = T))
tm_shape(regions2) +
tm_fill()
what is going on in terms of geometries? Behind the scene, both aggregate and summarize combine the geometries and dissolve the boundaries them using st_union
us_west = us_states[us_states$REGION == "West", ]
us_west_union = st_union(us_west)
texas = us_states[us_states$NAME == "Texas", ]
texas_union = st_union(us_west_union, texas)
## although coordinates are longitude/latitude, st_union assumes that they are planar
plot(st_geometry(us_west))
plot(texas_union)
plot(us_west_union)
st_cast behave differently on single simple feature geometry sfg objects, sfc, and simple feature objects
multipoint <- st_multipoint(matrix(c(1, 3, 5, 1, 3, 1), ncol = 2))
# st_cast can useful to transform the new object into linestring or polygon
linestring <- st_cast(multipoint, "LINESTRING")
polyg <- st_cast(multipoint, "POLYGON")
# reverse
multipoint_2 = st_cast(linestring, "MULTIPOINT")
multipoint_3 = st_cast(polyg, "MULTIPOINT")
all.equal(multipoint, multipoint_2, multipoint_3)
## [1] TRUE
multilinestring_list = list(matrix(c(1, 4, 5, 3), ncol = 2),
matrix(c(4, 4, 4, 1), ncol = 2),
matrix(c(2, 4, 2, 2), ncol = 2))
multilinestring = st_multilinestring((multilinestring_list))
multilinestring_sf = st_sf(geom = st_sfc(multilinestring))
multilinestring_sf
## Simple feature collection with 1 feature and 0 fields
## geometry type: MULTILINESTRING
## dimension: XY
## bbox: xmin: 1 ymin: 1 xmax: 4 ymax: 5
## epsg (SRID): NA
## proj4string: NA
## geom
## 1 MULTILINESTRING ((1 5, 4 3)...
linestring_sf2 <- st_cast(multilinestring_sf, "LINESTRING")
linestring_sf2
## Simple feature collection with 3 features and 0 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 1 ymin: 1 xmax: 4 ymax: 5
## epsg (SRID): NA
## proj4string: NA
## geom
## 1 LINESTRING (1 5, 4 3)
## 2 LINESTRING (4 4, 4 1)
## 3 LINESTRING (2 2, 4 2)
Geometric raster operations include the shift, flipping, mirrowing, scaling, rotation or wraping of images.
How to extract values from a raster overlaid by oter saptial objects. To retrive a spatial output, we can use almost the same subsetting syntax. The only difference is that we have to make clear that we would like to keep the matrix structure by setting the dop -parameter to FALSE
data("elev", package = "spData")
clip <- raster(xmn = 0.9, xmx = 1.8, ymn = -0.45, ymx = 0.45,
res = 0.3, vals = rep(1, 9))
elev[clip, drop = F]
## class : RasterLayer
## dimensions : 2, 1, 2 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : 1, 1.5, -0.5, 0.5 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 18, 24 (min, max)
When merge or performing map algebra on rasters, their resolution, project, origin and extent have to match. The same problems arises when we would like to merge satellite imagery from different sensors with different projections and resolutions. We can deal with such mismatches by aligning the rasters.
In the simplest case, two image only differ with regarding to their extent. Following code adds one row and two columns to each side of the raster while setting all new values to an elevation of 1000 meters
data(elev, package = "spData")
elev_2 <- extend(elev, c(1, 2), value = 1000)
map(list(elev, elev_2), extent)
## [[1]]
## class : Extent
## xmin : -1.5
## xmax : 1.5
## ymin : -1.5
## ymax : 1.5
##
## [[2]]
## class : Extent
## xmin : -2.5
## xmax : 2.5
## ymin : -2
## ymax : 2
Performing an algebra operation on two objects with differing extents in R, the raster packages returns the result for the intersection, and says so in a warnings.
elev_3 <- elev + elev_2
## Warning in elev + elev_2: Raster objects have different extents. Result for
## their intersection is returned
plot(elev_3)
identical(elev_3, raster::intersect(elev, elev_2))
## [1] FALSE
However, we can also align the extent of two rasters with extend. Instead of telling the function how may rows or columns should be added. We allow it to figure it out by using another raster object. Here, we extend the elev object to the extent of elev_2.
elev_4 <- extend(elev, elev_2)
the origin or a raster is the cell cornor closest to the coordinate (0,0). You can change raster orgin if two raster have different origins, since different origin make cells do not overlap completely which make map algebra impossible.
origin(elev_4) <- c(0.25, 0.25)
plot(elev_4)
plot(elev, add = T)
Different dataset can also differ with regard to their resolution. To match resolutions, one can either decrease or increase the resolution of one raster. Basically match resolution between different raster objects.
data("dem", package = "RQGIS")
dem_agg <- aggregate(dem, fact = 5, fun = mean)
plot(dem)
plot(dem_agg)
dem_disagg <- disaggregate(dem_agg, fact = 5, method = "bilinear")
plot(dem_disagg)
identical(dem_disagg, dem)
## [1] FALSE
The process of computing values of new pixel locations is also called resampling. In fact, the raster pacakge provides a resample fucntion. It lets you align several raster properties in one go, namely orign, extent and resolution. By default, it uses the bilinear -interpolation
dem_agg <- extend(dem_agg, 2)
dem_disagg_2 <- resample(dem_agg, dem)
Finally, in order to align many (possible hundres or thoughsand of) iamges stroed on disk, you could use the gdaltils::align_rasters function. However, you may also use raster with very large dataset.
This section focuses on interactions between raster and vector geographic data models. it includes four main techniques:
Many geographic data objects involve integrating data from many different soruces, such as remote sensing image (raster) and adminstrative boundaries. Often the extent of input raster datasets is larger thant the area of interest. In these case raster cropping and masking are useful for unifying the spatial extent of input data.
srtm = raster(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = st_read(system.file("vector/zion.gpkg", package = "spDataLarge"))
## Reading layer `zion' from data source `/Library/Frameworks/R.framework/Versions/3.6/Resources/library/spDataLarge/vector/zion.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 11 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 302903.1 ymin: 4112244 xmax: 334735.5 ymax: 4153087
## epsg (SRID): NA
## proj4string: +proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
map(list(srtm, zion), crs) # see that two file have different crs, so we need to reporject them
## [[1]]
## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
##
## [[2]]
## NULL
zion = st_transform(zion, projection(srtm))
crop() reduces the rectangular extent of the object passed to its first argument based on the exent of the object passed to its second argument.
srtm_cropped <- crop(srtm, zion)
Related to crop() is the raster function mask(), which sets values outside of bounds of the object passed to its second argument to NA. change the setting of mask yield different results. Setting updatevalue = 0, for example, will set all pixels outside the national park to 0. Setting inverse = TRUE will mask inside the bounds of the park pixels to 0
srtm_mask <- mask(srtm, zion)
srtm_inv_mask <- mask(srtm, zion,inverse = T)
plot(srtm_cropped)
plot(srtm_mask)
plot(srtm_inv_mask)
point extraction
data("zion_points", package = "spDataLarge")
zion_points$elevation = raster::extract(srtm, zion_points)
Raster extraction also works with lines selectors.
zion_transect = cbind(c(-113.2, -112.9), c(37.45, 37.2)) %>%
st_linestring() %>%
st_sfc(crs = projection(srtm)) %>%
st_sf()
transect = raster::extract(srtm, zion_transect, along = TRUE, cellnumbers = TRUE)
the along= TRUE and cellnumber = TRUE arguments to return cell IDs along the path. the result is a list contains a matrix of cell IDs in the first column and elevation values in the second.
The subsequent code chunk first converts this tricky matrix-in-a-list object into a simple data frame, returns the coordiantes associate with each extracted cell, and finds the associated distances along the transect.
transect_df <- purrr::map_dfr(transect, as.data.frame, .id = "ID")
transect_coords <- xyFromCell(srtm, transect_df$cell)
pair_dist = geosphere::distGeo(transect_coords)[-nrow(transect_coords)]
transect_df$dist = c(0, cumsum(pair_dist))
ggplot(transect_df) +
geom_line(aes(x = dist, y = srtm)) +
labs(x = "Distance (m)", y = "Elevation (m a.s.l.)",
title = "Elevation along the line") +
theme_ipsum(grid = "XY")
The final type of geographic vector object for raster extraction is polygons, like lines and buffers, polygons tend to return many raster valeus per polygon.
zion_srtm_values <- raster::extract(x = srtm, y = zion, df = T)
such results can be used to generate summary statistics for raster values per polygon, for example to characterize a single region or compare many regions.
zion_srtm_values %>%
group_by(ID) %>%
summarize(min = min(srtm),
mean = mean(srtm),
median = median(srtm),
max = max(srtm))
## # A tibble: 1 x 5
## ID min mean median max
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1122 1818. 1836 2661
data("nlcd", package = "spDataLarge")
zion_nlcd = raster::extract(nlcd, st_transform(zion, projection(nlcd)), df = TRUE, factors = TRUE)
zion_nlcd %>%
select(ID, levels) %>%
gather(key, value, -ID) %>%
group_by(ID, key, value) %>%
tally() %>%
spread(value, n, fill = 0)
## # A tibble: 1 x 9
## # Groups: ID, key [1]
## ID key Barren Cultivated Developed Forest Herbaceous Shrubland
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 leve… 98285 62 4205 298299 235 203701
## # … with 1 more variable: Wetlands <dbl>
Fast way to extracting raster cell values from a range of input geographic objects.
Parallelization : using many geographic vector select object by splitting them into groups and extracing cell values independently for each group (raster::clusterR() for details of this approach). WHen we use extract in polygons extraction we automatic use this function
velox pacakge, which provides a faster method for extracting raster data that fit in memory
Using R-GIS bridges: efficient calcuation of raster statistics from polygons can be found in the SAGA function saga:gridstatisticforpolygons
Rasterization is the conversion of vector data object into their representation in raster objects.
The rasterize() function contain two arguments x, vector object to be rasterized and y, a ‘template raster’ object defining the extent, resolution and CRS of the output. the geographic resolution of the input raster has a major impact on the results.
cycle_hire_osm_projected = st_transform(cycle_hire_osm, 27700)
raster_template = raster(extent(cycle_hire_osm_projected), resolution = 1000,
crs = st_crs(cycle_hire_osm_projected)$proj4string)
Rasterization is a very flexible operation: the results not only on the nature of the template raster, but also on the type of input vector and a variety of arguments taken by the rasterize function
First, we create a raster representing the presence or absence of cycle hire points.
ch_raster1 <- rasterize(cycle_hire_osm_projected, raster_template, field = 1)
plot(ch_raster1)
the fun argument speficies summary statistics used to convert multiple observation in close proximity into assocaite cells in the raster object. By default fun = “last” is used but other options such as fun = “count” can be used
ch_raster2 <- rasterize(cycle_hire_osm_projected, raster_template, field = 1, fun = "count")
plot(ch_raster2)
the new output, ch_raster2, shows the number of cycle hire points in each grid cell. The cycle hire location have different numbers of bicycles described by the capacity variable, rasing the question, what’s the capacity in each grid cell.
ch_raster3 <- rasterize(cycle_hire_osm_projected, raster_template,
field = "capacity", fun = sum)
plot(ch_raster3)
california = dplyr::filter(us_states, NAME == "California")
california_borders = st_cast(california, "MULTILINESTRING")
raster_template2 = raster(extent(california), resolution = 0.5,
crs = st_crs(california)$proj4string)
Line rasterization is demonstrated in the code below.
california_raster1 <- rasterize(california_borders, raster_template2)
plot(california_raster1)
Polygon rasterization, by contrast, selects only cells whose centroids are inside polygon
library(microbenchmark)
california_raster2 <- rasterize(california, raster_template2)
california_raster3 <- fasterize(california, raster_template2)
microbenchmark(rasterize(california, raster_template2), fasterize(california, raster_template2))
## Unit: microseconds
## expr min lq mean
## rasterize(california, raster_template2) 60709.583 63782.450 68666.7704
## fasterize(california, raster_template2) 35.227 45.104 70.9925
## median uq max neval
## 66941.7635 70040.5965 84465.180 100
## 87.4395 92.0715 109.407 100
plot(california_raster2)
It involves converting spatially continous raster data into saptially discrete vector data such as points, lines or polygons
Keep in mind. In R, vectorization refers to the possibliiyt of replacing for- loops and alike by doing things like 1:10 / 2
elev_point <- rasterToPoints(elev, spatial = T)
plot(elev)
plot(elev_point)
Another common type of spatial vectorization is the creation of contour lines representing lines of continuous height or temperatures for example
data(dem, package = "RQGIS")
cl = rasterToContour(dem)
plot(dem, axes = FALSE)
plot(cl, add = TRUE)
Contours can also be added to existing plot with function such as contour an contourplot or tm_iso
# create hillshade
hs <- hillShade(slope = terrain(dem, "slope"), aspect = terrain(dem, "aspect"))
plot(hs, col = gray(0:100 / 100), legend = F)
# overlay with DEM
plot(dem, col = terrain.colors(25), alpha = 0.5, legend = F, add = T)
# add contour lines
contour(dem, col = "white", add = T)
The final type of raster vectoriization involves conversion of raster to polygons.
This is illustrate by converting the grain object into polygons and subsequently dissolving borders between polygons with the same attributes. Attributes in this case are stored in a column called layer
grain_poly <- rasterToPolygons(grain) %>%
st_as_sf()
grain_poly2 <- grain_poly %>%
group_by(layer) %>%
summarize()
# or just rasterToPolygons(grain, dissolve = T) %>% st_as_sf()
library(RQGIS)
## Loading required package: reticulate
data(random_points)
data(ndvi)
ch = st_combine(random_points) %>%
st_convex_hull()
nz_geo <- st_geometry(nz)
map(c(100, 1000, 10000, 100000), ~ plot(st_simplify(nz_geo, dTolerance = .x)))
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
plot(st_simplify(nz_geo, dTolerance = 100000, preserveTopology = T))
Use st_simplify function, it start to bread down at 100,000.
Answer: problems: st_simplify returns polygon and multipolygon resuts, affecting plotting. Cast into a simple geometry type to resolve this
nz_simple_multipoly <- st_simplify(st_geometry(nz), dTolerance = 10000) %>%
st_sfc() %>%
st_cast("MULTIPOLYGON")
plot(nz_simple_multipoly)
map(c(0.5, 0.005, 0.0005, 0.00005), ~ plot(ms_simplify(nz_geo, keep = .x)))
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
Use ms_simplify, nz start to break down at 0.05 level
nz_simply1 <- st_simplify(nz_geo, dTolerance = 100)
nz_simply2 <- ms_simplify(nz_geo, keep = 0.5)
map(list(nz_geo, nz_simply1, nz_simply2), class)
## [[1]]
## [1] "sfc_MULTIPOLYGON" "sfc"
##
## [[2]]
## [1] "sfc_GEOMETRY" "sfc"
##
## [[3]]
## [1] "sfc_MULTIPOLYGON" "sfc"
After simplify, result from st_simplify change to sfc_GEOMETRY instead of sfc_multipolygon.
canterbury <- nz %>% filter(Name == "Canterbury")
nz_height_buffer <- st_buffer(nz_height, dist = 100)
plot(st_geometry(nz))
plot(st_geometry(nz_height_buffer), add = T)
canterbury_within <- st_intersection(canterbury, nz_height_buffer)
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries
nrow(canterbury_within)
## [1] 75
There are 75 out of 101 highest points inside New Zealand within 100 km of canterbury.
nz_centroid <- st_centroid(st_union(nz))
canterbury_centroid <- st_centroid(canterbury)
## Warning in st_centroid.sf(canterbury): st_centroid assumes attributes are
## constant over geometries of x
st_distance(nz_centroid, canterbury_centroid)
## Units: [m]
## [,1]
## [1,] 234192.6
234192.6 meters, is the distance between center of nz and center of canterbury
world_sfc <- st_geometry(world)
world_southup <- world_sfc * matrix(c(-1,-1), ncol = 2)
plot(world_southup)
china_sfc <- world %>%
filter(name_long %in% c("China", "Taiwan")) %>%
st_union() %>%
st_geometry()
china_southup <- china_sfc * matrix(c(-1,-1), ncol = 2)
plot(china_southup)
plot(box)
plot(x, add = TRUE)
plot(y, add = TRUE)
plot(p, add = TRUE)
text(x = c(-0.5, 1.5), y = 1, labels = c("x", "y"))
p[b]
## Geometry set for 8 features
## geometry type: POINT
## dimension: XY
## bbox: xmin: -0.882033 ymin: 0.7773639 xmax: 1.772728 ymax: 1.882
## epsg (SRID): NA
## proj4string: NA
## First 5 geometries:
## POINT (1.772728 1.348663)
## POINT (-0.1341215 0.8641556)
## POINT (1.310264 0.9987838)
## POINT (1.318306 0.7773639)
## POINT (-0.882033 0.7907506)
st_intersection(b, p)
## Geometry set for 9 features
## geometry type: POINT
## dimension: XY
## bbox: xmin: -0.882033 ymin: 0.7773639 xmax: 1.772728 ymax: 1.882
## epsg (SRID): NA
## proj4string: NA
## First 5 geometries:
## POINT (1.772728 1.348663)
## POINT (-0.1341215 0.8641556)
## POINT (1.310264 0.9987838)
## POINT (1.318306 0.7773639)
## POINT (-0.882033 0.7907506)
us_states %>%
mutate(state_border = st_length(geometry)) %>%
arrange(desc(state_border))
## Simple feature collection with 49 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -124.7042 ymin: 24.55868 xmax: -66.9824 ymax: 49.38436
## epsg (SRID): 4269
## proj4string: +proj=longlat +datum=NAD83 +no_defs
## First 10 features:
## GEOID NAME REGION AREA total_pop_10 total_pop_15
## 1 48 Texas South 687714.3 [km^2] 24311891 26538614
## 2 06 California West 409747.1 [km^2] 36637290 38421464
## 3 26 Michigan Midwest 151119.0 [km^2] 9952687 9900571
## 4 12 Florida South 151052.0 [km^2] 18511620 19645772
## 5 30 Montana West 380829.2 [km^2] 973739 1014699
## 6 27 Minnesota Midwest 218566.5 [km^2] 5241914 5419171
## 7 16 Idaho West 216512.7 [km^2] 1526797 1616547
## 8 51 Virginia South 105404.7 [km^2] 7841754 8256630
## 9 35 New Mexico West 314886.1 [km^2] 2013122 2084117
## 10 53 Washington West 175436.0 [km^2] 6561297 6985464
## state_border geometry
## 1 4957325 [m] MULTIPOLYGON (((-103.0024 3...
## 2 3793650 [m] MULTIPOLYGON (((-118.6034 3...
## 3 3579186 [m] MULTIPOLYGON (((-85.63002 4...
## 4 2958238 [m] MULTIPOLYGON (((-81.81169 2...
## 5 2827680 [m] MULTIPOLYGON (((-116.0492 4...
## 6 2567879 [m] MULTIPOLYGON (((-97.22904 4...
## 7 2567784 [m] MULTIPOLYGON (((-116.916 45...
## 8 2406526 [m] MULTIPOLYGON (((-75.66971 3...
## 9 2378365 [m] MULTIPOLYGON (((-109.0452 3...
## 10 2346224 [m] MULTIPOLYGON (((-122.7699 4...
Texas have the highest border line.
plot(ndvi)
plot(st_geometry(random_points), add = T)
plot(ch, add = T)
# crop, reduce the rectangular extent of the object passed to its first argument based on the extent of the object passed to its second argument
ndvi_ran_crop <- crop(ndvi, as(random_points, "Spatial"))
spplot(ndvi_ran_crop)
ndvi_ch_crop <- crop(ndvi, as(ch, "Spatial"))
spplot(ndvi_ran_crop)
# mask compare to crop, preserve the extent of the first object, but sets values outside of the object passed to its second argument to NA
ndvi_ran_mask <- mask(ndvi, as(random_points, "Spatial"))
spplot(ndvi_ran_mask)
ndvi_ch_mask <- mask(ndvi, as(ch, "Spatial"))
spplot(ndvi_ch_mask)
random_buffer <- st_buffer(random_points, dist = 90)
plot(ndvi)
plot(st_geometry(random_buffer), add = T)
plot(ch, add = T)
random_points$ndvi1 <- raster::extract(ndvi, as(random_points, "Spatial"))
random_points$ndvi2 <- raster::extract(ndvi, as(random_points, "Spatial"), buffer = 90, fun = mean, na.rm = T)
ggplot(random_points) +
geom_point(aes(x = ndvi1, y = ndvi2)) +
theme_ipsum_es()
Extracting values by buffers is more suitable than by points alone, for example when we expect errors in single raster cells. Average values can decreases an effect of errorness values in some raster vells
nz_3100 <- nz_height %>%
filter(elevation > 3100)
nz_3100_template <- raster(extent(nz_3100), resolution = 3000,
crs = st_crs(nz_3100)$proj4string)
plot(st_geometry(nz_3100), graticule = st_graticule(nz_3100, datum = 2193), axes = T)
high_raster1 <- rasterize(nz_3100, nz_3100_template, field = 1, fun = "count")
spplot(high_raster1)
high_raster2 <- rasterize(nz_3100, nz_3100_template, field = "elevation", fun = "max")
spplot(high_raster2)
high_raster3 <- raster::aggregate(high_raster1, fact = 2, fun = sum)
spplot(high_raster3)
high_raster4 <- resample(high_raster3, high_raster1)
spplot(high_raster4)
Advange: low memory use advantage: faster processing advantage: good for vis in some cases Disadvange: remove geographic detail
grain_poly <- rasterToPolygons(grain) %>%
st_as_sf()
clay <- grain_poly %>%
filter(layer == 1)
plot(clay)
Advange: can be used to subset other vector object, can do affine transformation and use sf/dplyr verbs Disadvange: better consistency, fast processing on some operation.